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Abstract 

The construction of a statistical model for eigenfunctions of the Ising model in transverse and longitudinal fields is discussed 
in detail for the chaotic case. When the number of spins is large, each wave function coefficient has the Gaussian distribution 
with zero mean and the variance calculated from the first two moments of the Hamiltonian. The main part of the paper 
is devoted to the discussion of different corrections to the asymptotic result. One type of corrections is related with higher 
order moments of the Hamiltonian and can be taken into account by Gibbs-like formulae. Another corrections are due to 
symmetry contributions which manifest as different numbers of non-zero real and complex coefficients. Statistical model with 
these corrections included agrees well with numerical calculations of wave function moments. 
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I. INTRODUCTION 


Direct numerical calculations of non-integrable quantum many-body problems become exponentially difficult for 
large number of particles. The developing of approximate statistical description of such models has been and remain 
thus of great importance. 

For macroscopic bodies the number of interacting particles, N, is so large that only the thermodynamic limit 
N ^ oo is meaningful. In this case the behaviour of a small sub-system is described by its reduced density matrix 
constructed from knowledge of eigenvalues, E^, and eigenfunctions, ipin), when the interaction with outside particles 
is switched off. According to general properties of quantum thermodynamics (see e.g. [T]) the mean value of an 
observable A of the sub-system is calculated as follows 

1™;^/ (^'(t)|i|«'(t))dt = ^^e-'^'®"(^/>(n)|i|i/'(n)), Z = (1) 

T-i-OO T In Z ^ ^ 


and the only information of the existence of a larger system is the inverse temperature /? which is determined by the 
mean energy of the initial state. 

This text-book approach reveals so successful that there is no doubt (at least in physical literature) that it is the 
correct way of calculation of the thermodynamic limit in generic systems. The existence of powerful perturbation 
expansion around equilibrium states with zero and non-zero temperature [2] , and non-equilibrium states [3] contributes 
considerably to the success of quantum thermodynamics. The important particularity of these constructions is that 
in any steps effectively infinite number of particles never appears explicitly. The formalism is built in such a way that 
only finite quantities like temperature and chemical potential give information of the existence of outside world. 

With the ever increasing power of modern computers, it has become possible to perform full numerical calculations 
for quantum systems containing few tens of particles. This in turn enables meaningful formulations of questions 
about the foundations of quantum thermodynamics, which is a work in progress in modern physics: what systems are 
thermalized, what are corrections to thermodynamic limit, etc. 

Let Ea and 'I'(q;) be exact eigenvalues and eigenfunctions of the full N particles Hamiltonian. For finite N, instead 
of canonical Gibbs measure 0. one has to use the microcanonical average 


lim — 
T-»-oo T 


(vI/(t)|A|vI/(t))dt = 




En.AEn 


(«'(a)|A|^'(a)) 


\Eo — Ea\<AEo 


( 2 ) 


where Eq is the mean energy of the initial state, AEq is a energy window assumed to be small with respect to Eq but 
much larger than the distance between nearest levels 


l/p(A'o) AEq ^ Eq (3) 

where p{E) is the mean spectral density at energy E. NeqAEo is the number of states in this window. (We consider 
only the so-called chaotic systems where there is no integral of motion except classical ones.) 

The important difference between Eq. Q and Eq. ([^ is that in the former ’tp{n) is eigenfunctions of a sub-system 
but in the later 'l'(a) are eigenfunctions of the full Hamiltonian. 

In such approach the central object of investigation is A^-particles eigenfunction ^'(a). There exist different scenarios 
about the emergence of thermodynamic behaviour (see e.g. [4]). One of the most accepted is the so-called eigenstate 
thermalization conjecture (ETH) [5]-[7j which states that relation 0 is valid for (almost) all individual eigenstates. 
Similar conjecture is known in low-dimensional quantum chaos (see e.g. i) where it is explicitly stated that for a 
smooth classical observable A(p, q) almost all diagonal matrix elements of the corresponding quantum counterpart 
tend to the microcanonical semiclassical limit when h ^ 0 

{'^{a)\A\'i>{a)) J 6{E - H{p,q))A{p,q)dpdq, Z = J 6(E - H{p,q))dpdq. (4) 

The smallness of off-diagonal matrix elements (4 'q,|A|4'^) ^ (4'q,|A|4'q,), when ^ a also appears naturally in 
quantum chaos (cf. e.g. [5]) as only periodic orbits give the contribution to this quantity. 

These and many other arguments strongly suggest that for systems where quantum thermodynamics can be applied, 
the majority of eigenfunctions are universal depending only on a few parameters. It has been even proved that in (a 
certain sense) typical functions in quantum mechanics lead automatically to canonical averaging |13| . |I4j . The well 
known example of such phenomenon is given by Berry’s conjecture |in)-|I2j which postulates that wave functions of 
low-dimensional chaotic systems are Gaussian random functions with variance determined from semiclassical micro- 
canonical average as in Eq. 0. Similar construction has been applied for various problems in nuclear and atomic 
physics, many-body models etc, 0, iiij-[5n]. 
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The purpose of this paper is to construct and carefully check a statistical model for eigenfunctions of the quantum 
Ising model in transverse and longitudinal fields determined by the Hamiltonian 


N 


N 


N 




(5) 


n—1 


n—1 


n—1 


where are usual Pauli matrices, and parameters A and a fix values of transverse and longitudinal fields. 

This model is a prototypical example of one-dimensional spin chains with nearest-neighbour interactions. When 
a = 0 it reduces to the well known quantum Ising model in transverse field which is integrable by the Jordan-Wigner 
transformation |31] . [32j and for A = 1 it becomes critical and serves as paradigmatic model of quantum critical 
phenomenon |34j . For non-zero a (and A 7 ^ 0) the model is considered as non-integrable but at critical value of 
transverse field, A = 1, and specially fine-tuned weak longitudinal field, a —> 0, it is integrable but not conformal [33]. 

In |35| it was shown that the ground state wave function of the Ising model (as well as of practically all one¬ 
dimensional spin-chain models) is multifractal in the initial spin basis. So it does not feet the standard thermodynamic 
scheme. For high excited states the situation is different. The spectral density for this model has been discussed in 
[36j . For finite number of spins, iV, in the bulk of the spectrum when E ^ -s/iV there exist two different regimes. 
When the both fields A and a are of the order of 1 (i.e. of order of the hopping term), the spectral density is well 
approximated by a simple Gaussian function whose parameters are calculated from the knowledge of the first moments 
of the Hamiltonian. If A is small or large the spectral density of the Ising model in two fields for large but finite N 
has many peaks well described by a sum of Gaussian functions calculated directly from the Hamiltonian without full 
diagonalization. 

In parallel with two regimes of spectral density of the Ising model there are two regimes for eigenfunctions in 
the bulk of the spectrum. Here we consider in detail the case when all coupling constants are of the same order. 
Multi-peaks case will be discussed elsewhere. 

The plan of the paper is the following. The construction of chaotic wave functions for the Ising model is discussed 
in Section [TT| To check the accuracy of such statistical model it is convenient to calculate moments of wave functions 
and compare them with statistical approximations. In Section HI it is done for the participation ratio for the Ising 
model in two fields. The main part of this Section is devoted to the calculations of different corrections to the lowest 
order approximation. It is demonstrated that corrected statistical model of eigenfunctions agrees very well with the 
results of direct numerical calculations. The first moments of the full Ising model Hamiltonian are determined in 
Appendix 


II. STATISTICAL DESCRIPTION OF WAVE FUNCTIONS 

The construction of chaotic wave functions for different problems follow approximately the same steps as in m- 

An eigenfunction of the Hamiltonian ([^ with energy Ea by definition obeys 

H|^„) = (6) 

and can be represented as a finite series of basis set functions, \n) 

|'ko)=E^"(«)l«)- (7) 

n 

For N spins-i the full dimension of the Hilbert space is V = 2^. Here a state is indicated by symbol ft = 
(ni, 712,..., u-at) with the convention that Uk = 0 corresponds to spin at site k down and rifc = 1 shows that site 
k is occupied by spin up. 

All information of wave function is contained in coefficients Cfi{a). Let us consider the collection of coefficients 
Cfj(a) with fixed symbol ft and with energies in a small interval I = [E < E^ < E + AE], The energy window AE is 
assumed to be small with respect to E but much larger than the distance between nearest levels as in ([^. 

The principal assumption is that for chaotic systems coefficients C's(a) in such intervals are so irregular functions 
of eigen-energies E^ that their exact values are irrelevant for most purposes and one has to developed statistical 
description of chaotic wave functions. For spin-chains models without random parameters as the one given by Eq. (§ 
we are not aware of strict definition of chaoticity. The usual lore is that non-integrability of systems with large number 
of degrees of freedom is (almost) synonym of chaoticity. For our purposes such notion is sufficient. 

As an example, we present in Fig. one particular coefficient for all eigen-energies of the Ising model with a = 1 
and A = 1 for 17 spins with periodic boundary conditions obtained by direct numerical diagonalization. This type of 
pictures strongly suggests the following conjecture. 
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FIG. 1. (a) Coefficient corresponding to symbol | 00001000110111001) for the Ising model with 0 = 1 and A = 1 for = 17 
spins in the sector with zero translational momentnm. Red vertical lines indicate energy windows used in the calculation of 
distributions of this coefficient in Fig. [^). Each window contains approximately 2000 levels, (b) Distribution of the coefficient 
as in Fig. &) in 4 energy windows indicated in that figure (red histograms). Blue solid lines are the best Gaussian fits to these 
distributions. Inserts: The difference between numerical histograms and the Gaussian fits. 


Conjecture 1: The irregular behaviour of coefficients Cfj(a) in a small energy window can be mimic by the 
assumption that they are random functions of energies with a certain distribution. 

For strongly chaotic systems there are two widely used a-priori assumptions. If coefficients are real, they should well 
be described by the Gaussian distribution. It means that at fixed symbol n, the distribution density of coefficients is 


P{Cn = x) 




( 8 ) 


If coefficients are complex, one assumes that they have Gaussian distribution for the both real and imaginary parts 
with zero means and the same variance 


P(ReCf; = x.luvCn = y) 


1 f x'^ +y'^\ 
""P [ E| j ■ 


In the both expressions indicates the mean value of modulus square of coefficient Cfi 


(9) 


S| = {\Cn?) (10) 

For illustration, in Fig. we compare distribution of the same coefficient as in Fig. &) with the best Gaussian fit 
and a good agreement is clearly seen. 

Accepting the above conjecture means, in particular, that moments of variables Cfj averaged over a small en¬ 
ergy interval I should be well approximated (at least for large N) hy the moments of the corresponding Gaussian 
distributions 




E 


E^ei 


where for complex coefficients TZg = fQj. j-gg,! ones TZg = TZg^^'^ with 


^(complex) ^ ^(real) ^ 


2«r {q 




( 11 ) 


( 12 ) 


In particular, _ 2 and = 3. 
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Under the validity of the conjecture all average moments of Cfi{a) (with fixed symbol n) are determined by one 
quantity, the variance S| = which depends on symbol n and of the center of energy window E. 

To calculate this variance, it is convenient to consider the so-called strength function (or the local density of states) 


Pft{E) =Y,\Cuia)\^SiE - E^). 

a 


(13) 


Due to pseudo-random character of Cfi{a) and E^ the strength function itself can be considered as pseudo-random 
or random function of energy E. The mean value of a certain function f{E) is defined as in © 


iPn{E)f{E)) 


1 

J^E.AE 


^ |C„-(a)|V(i?a) 

\E-Ea,\<AE 


(14) 


where the width of the energy window, AE, obeys 

The knowledge of the strength function permits to find \Cft{a)\'^ averaged over a small energy window (which can 
be used as the definition of the variance) 


,2 {Pnm 

PiE) 


(15) 


where p{E) is the average density of states with energy E i.e. the number of states in an interval E < Ea < E + AE. 
The advantage of the strength function is that for this quantity there exist exact sum rules m 


J Pft{E)E>^dE = {E^} = {n\n^\n) 


(16) 


It means that these moments can be calculated directly from the Hamiltonian without solving the full problem. 

Of course, the exact calculation of the strength function is equivalent to the full solution of the problem. To obtain 
a simple approximation one has to assume that the strength function can be well approximated from the knowledge 
of its first moments. One can argue that for iV-body systems with short-range interactions the first two moments 
reproduce well higher order moments when N —>■ oo. It leads to the second conjecture: 

Conjecture 2: The functional dependence of Pn{E) on E is well approximated when —>■ oo by a Gaussian 


PfiiE) = 


Sttct? 


: exp - 


(E - Enf 




where Efj and cr- are the first two moments of the Hamiltonian 

En = {fi\H\n), = {n\{n - Enf\n) 

Under this conjecture the variance of coefficient distribution given by Eq. (|15[) is 


1 


p{E) J2-Ka. 


: exp - 


{E - Eaf 




(17) 


(18) 


(19) 


with the mean spectral density calculated by the formula (to ensure the normalization \Cfi{a)\‘^ = 1) 




27r CT? 


: exp - 


{E - Enf 




Assuming that coefficients have Gaussian distribution the values of 2q*‘'^ moments of the wave function are 


( 20 ) 


Af,s(y:ic„(o)p") = ^y: 




exp -q 


(E - Eaf 


2al 


( 21 ) 


The above conjectures are basis ingredients of construction of statistical models for chaotic wave functions in different 
problems H], [I5]-[3D]. Though they were not proved in full generality (but see [37], [33]), they are simple enough to 
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be check for a given particular system. One can even reverse the arguments and say that wave function of A^-body 
model with short range interactions is called chaotic iff it obeys the above conjectures. 

Another line of reasoning may be the maximum-entropy principle (see e.g. |89j 1 according to which ’’the best” 
statistical distribution {Pft{E) in our case) fulfilled certain restrictions (as in & is given by the one which maximise 
the entropy 


S = - J ME) \nPfi{E) dE 


( 22 ) 


Assuming e.g. that a finite number of first moments in (16) are specified, the maximum-entropy principle predicts 


that the corresponding probability density takes the Gibbsian form 

3- E, 


Pn{E) = 


Zip) 


Z{pi)= J 


(23) 


with Lagrangian multipliers, calculated from the partition function using the knowledge of the first 

moments 


{n\W\n) =- In Z(fl) 


(24) 


Of course, when only two moments are taken into account one gets Hzl)- 

Though in all considered cases expression (23) gives the best results, it requires numerical calculations of Lagrangian 


multipliers which complicates the separation of different sources of corrections. To increase the accuracy of approxi¬ 
mation, we shall include the third and the forth moments assuming that that they are much smaller than the first two 
moments by using Gram-Charlier A series in Hermite polynomials. It leads to the following modification of Pn(E) 


Pu{E) = 


1 




: exp 


TTCr-^ 


{E - Er,f 

2 a 2 




/ ^ \ TT [ ^ \ 

I ct ) 


E — En 


3!cr^ 


(25) 


Here and are the third and the forth cumulants for the Ising model given by Eqs. (A18) and (A24| from 
Appendix [ a{ E[s{x) = — 3x and Hi{x) = x'^ — dx^ -\- 3. 


III. PARTICIPATION RATIO FOR THE ISING MODEL IN TWO FIELDS 


For the Hamiltonian © simple calculations give 

En = \{N -2n), al = a'^ = N{l + a^) 

where n is the number of spins up. Each state with fixed n is C]^ degenerated. 
These formulae mean that 

P 1 ( (E-E^)^ 

Pn{E) = ^—_ exp ( - 


\/2Tra'^ 


2 a 2 


(26) 


(27) 


Notice that we wrote n and not n to stress that in this approximation all quantities depend only on number of spins 
up. 

The full density normalized to 1 is 


p{E) = 2 


N 

-N 


y-in p 

n =0 


{E) 


(28) 


In Fig. ([^ this equation is compared with the spectral density of the Ising model ([^ with a = 1 and A = 1 for A = 17 
spins calculated numerically. The agreement is good but small differences are visible. When higher moments taken 
into account as in (25) the result is practically indistinguishable from numerics. 


To approximate moments of wave functions accurately one has to know not only the first moments of wave function 
coefficients but also are they real or complex. Usually, the answer is simple. If the Hamiltonian is real symmetric, 
wave functions (and their coefficients) are real. If the Hamiltonian is complex Hermitian, coefficients are also complex. 
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FIG. 2. Spectral density for the Ising model q = 1 and A = 1 for A'’ = 17 spins. Blue histogram indicates the numerical results. 
Black dashed line is the Gaussian approximation without corrections (281 with (271. Red solid line the corrected approximation 
(281 with (251. 


For the considered Ising model in two fields ([^ the situation is more tricky. The point is that this model with periodic 
boundary conditions has the translational invariance 

T: n —?■ n + 1, n = l,2,...,N (29) 

so if is a wave function with energy E,^ = J2nj=±i Cni,n2,...,ni^\ni,n2, .. ■, un) then TT = J2nj=±i C'ni.n 2 .....nN ni,..., un¬ 

is also an eigenfunction with the same energy which imposes certain relation on coefficients. 

As = 1, its eigenvalues are and eigenfunctions can be classified by the translational momentum k = 


f (30) 

Numerical calculations below are performed in the basis with fixed translational momentum k. 

To construct explicitly such basis one fixes one element, |n), n = (ni, 712 ,..., uat) and constructs all its translations 
with corresponding phases. For prime N one can choose 


Here T-^ ln) indicates the shift on j elements of the initial sequence 

T-^ Ini, 712,..., 71 jv) = |777v-i+i, • ■ ■, nN, ui, 712, • ■ ■, nN-j) 


(32) 


After performing one step one takes arbitrary element which is not in the sum (311 and repeats the construction till 
all elements are exhausted. 


For composite N there exist elements with primitive periods t equals to a divisor of N. Their number is 


r(t) = i^2*/V(d) 


(33) 


d\t 


where is the Mdbius function, fJ-{n) = (—I)"*’ if n is square-free number with Up prime divisors, and n{n) = 0 if 
n is divisible on a square of a prime. 

The total dimension of the basis with fc = 0 is 

A4ot = 5] rin) = 1 ^ (34) 

n|Ar d\N 


where <(>(71) is Euler’s totient function equals the number of positive integers smaller than n and co-prime with it. For 
non-zero k the contribution of basis states with primitive period t < N, (t\N), is non-zero only when k = 0 mod N/t. 
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The simplest approximation 


oAf 

■M;ot ~ (35) 

is sufficient in many cases when N is large. The total number of states with fixed momentum and fixed number of 
spins up t'tot(^) can also be approximated in the similar manner 

i^totin) ~ (36) 

After rewriting the Hamiltonian in the basis of states with fixed translational momentum, |nfe), it becomes complex 
and eigenfunctions 


^'fc — 'y ( C'nfc 
fik 


(37) 


are, in general, also complex. 

Nevertheless, one cannot conclude that in the chaotic regime eigenstates and eigenfunctions are distributed as for 
Gaussian Unitary ensemble of random matrices (GUE). The reason is that for the model considered there exists 
another discrete symmetry, namely geometric inversion 

S:n^N-n + l, n=l,2,...,N (38) 

Hamiltonian (§ is invariant under this transformation. Therefore if is an eigenfunction S'l' is also an eigenfunction 
with the same energy. 

By definition the initial basis |n) transforms under S as follows 

S\ni,n2,...,nN) = \nN,nN-i, ■ ■ ■ ,ni) (39) 

Functions with zero translational momentum transform into themselves under this inversion 

(40) 

where e = ±1 is the parity under the inversion. But functions with k ^ Q are transformed up to a phase into functions 
with opposite momentum 

-54'fc = 4'Ar-fe (41) 

States |n)fe with hxed momentum k (cf. Eq. ©) split into two groups under the inversion. The Hrst group includes 
all states which under this inversion transform into themselves but with opposite momentum (once more up to a 
certain non-dynamical phase which can be included in the definition of basis states) 

S\n)k = \n)N-k (42) 

In other words, the inversion of one element of this group is equivalent to a certain shift of this element. Such elements 

are called invariant under the inversion. 

Elements of the second group are organised in different pairs transformed one to another by the inversion 

S\n)k = \n')N-k, fi'^n (43) 

We call them non-invariant elements. 

The above arguments show that eigenfunctions with k ^ 0, N/2 can be written in the form 

'1''.= E (clr‘°U)*+cU'"’'i«')‘)+ E M 

non—invariant invariant 


where n and n' are pairs of elements connected by the inversion. Taking real and imaginary parts of this expression 
demonstrates that all coefficients in the expansion can be chosen real which explains that spectral statistics of this 
model is well described by GOE statistics (either give the reference or a picture of left-right asymmetry for non-zero 
momentum!) though the Hamiltonian itself is complex (cf. |12]). 


For k = 0 eigenfunctions have a particular parity under the inversion and all coefficients can be chosen real but 
invariant elements are identically zero for states with negative parity 




= E 


ct"^\n)o + eCt"^\n' 


non—invariant 


)o) 


+ ^ 
^ 2 


E 

invariant 




(45) 


with e = ±1. 

It is plain that the majority of elements at large N are non-invariant. Nevertheless, corrections from invariant 
elements, though small, are noticeable for N accessible in numerical calculations. The number of invariant elements 
is especially simple to find for prime N when all elements have the same primitive period. As N = 2K -|- 1 is odd, 
all invariant elements in the periodic cycle are appeared in pairs under the inversion except one element which has 
to be invariant under the inversion. Such element can be constructed by choosing arbitrary the first K elements and 
reflecting them into another part of the element. The number of possibilities is 2^ and one has 2 chooses for the 
central element, therefore the total number of invariant elements is 

Afinv = 2^+1 = (46) 

where [x] is the integer part of x. The number of invariant elements with hxed total number of spins is 

(47) 


Similar expression may be obtained for even N. 

According to the above statistical conjectures, coefficients with k ^ 0 are distributed as complex Gaussian 


variables with zero mean and variance calculated from PniE) as in Eq. (15) 




PftiE) 
P{E) ■ 


N 


p{E) = ^ Vtot{n)Pn{E) 


(48) 


n=0 


with PniE) as in Eqs. (27) or (25). 


Similarly, coefficients of invariant terms, have real Gaussian distribution with the same variance. Therefore, 

when moments are calculated, one has to take into account corrections from different types of distributions (cf. 
Eqs. ([Hi). For large N one can neglect the existence of invariant elements and consider that for non-zero k functions 
are complex and for zero momentum functions are real. 

As a typical example we consider the participation ratio determined as usual 


Pr{E^) = (^\CniE^)\^'^ 


(49) 


The simple Gaussian approximation of this quantity takes the form 


Pr{E) « A4c 


^2 {‘2-^En=oC^PniEy 


(50) 


where PniE) is given by Eq. (27). Here we use approximation (36) for the number of total terms with fixed n. 
A/tot = ^nj^totin) = 2^/N, 77.2 = 2 for states with non-zero momentum and 77-2 = 3 for zero momentum state. 

In Figs. 1^) and j^) we present the numerical results for this quantity for the Ising model with a = A = 1 and 
N = n spins for states with translational momentum k = 2 and fc = 0 respectively. Notice that left-hand side of 
the participation ratio (with E < —10) fluctuate strongly and the using statistical description at energies close to the 
ground state is questionable. In the same figures the above Gaussian approximation is also indicated. In general the 
agreement is good but there exit noticeable differences with numerical calculations. 

There exist two different types of corrections. The first is related with higher moments of Hamiltonian and they 
can be taken into account by modifications indicated in Eq. (25). This formula differs from the Gibbs-like formulae 


(23). The latter usually give better results, but as they require numerical calculations to find necessary Lagrangian 


multiplies, their use leads to non-transparent expressions where relative importance of different terms is hidden. To 
get clear separation of different contributions we prefer to use series expansion as in Eq. (25) though they lead to big 


errors at large arguments. As in all cases the advocated statistical approach can, strictly speaking, be applied only in 
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(a) 



(b) 


FIG. 3. (a) Participation ratio of the Ising model with a = 1 and A = 1 for = 17 spins in sector with translational moment 
equals 2. Black crosses are result of numerical calculations. Blue dashed line is the pure Gaussian approximation without 
corrections (|50[l with 77.2 = 2. Red line is the corrected approximation, (b) The same as in (a) but with zero translational 


moment. Blue dashed line is the pure Gaussian approximation without corrections (501 but with 77.2 = 3. 


the bulk of the spectrum it is not an important restriction. Another corrections are related with the fact that in the 
expansion of eigen functions there exits basis sets which are invariant under inversion (cf. Eqs.(44) and ( |45[ ))which 
has properties different from non-invariant terms. 

For the moments of wave function coefficients with non-zero momentum one gets the following formula 


^ |4's(A)| 29 \ = (T^complex j^real _ 


(Elo^tot(n)P„(F))‘ 


where 


^,{E) = 


X/n ^inv(^) PnjE) ^ J^inv 
J2n ^tot (n) P^ (E) Aftot 


(51) 


(52) 


The approximative expression of « <5 = A/inv/Atot corresponds to the uniform approximation where the summation 
over invariant set is proportional to the summation over all elements. When the invariant elements are taken into 
account in Eq. (50) for the participation ratio one has instead of TZ 2 = 2 to use TZ 2 = 2 + S ioi states with non-zero 


translational momentum. 

For fc = 0 all eigenfunctions are either symmetric or anti-symmetric with respect to inversion. The full dimensions 
of symmetric and anti-symmetric sub-spaces (i.e. the number of independent real coefficients) are (cf. mi) 


iV± = i(l±(5)A4ot 


(53) 


Eigenvalues corresponding to the both sub-spaces are statistically independent, so the local densities of these eigen¬ 
values are proportional to the above numbers, 


P±{E) = ]^{l±5)p{E) 


(54) 


Correspondingly, the strength function (13) has contributions from the both symmetric and anti-symmetric states. 


According to the above conjecture each term in Eq. (45) has real Gaussian distribution with the same variance. But 


this variance are depended on the number of independent components. In the uniform approximation variances in 
different sub-spaces are 


For positive parity states the variance of equals S| (-I-) but variance of non-invariant coefficients, C- 


S|o(±) = 


_ 

l=b5 


(55) 


(n-in) 


no ' 


IS, 


according to our definition (45), equal to iE| (-I-). Calculating the 2q^^ moments of eigen-functions in the initial spin 
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FIG. 4. Inverse moments of the Ising model with a = 1 and A = 1 for = 17 spins with zero translational momentum for 
different values of q. Black crosses are results of numerical calculations. Red lines are the corrected Gaussian approximation. 


basis gives 


M(+) 


1 -( 5 + 29(5 
( 1 + 5)9 


Mg, 


Mi") = 


(1-5)' 




(56) 


As we do not separate states with definite parity, the mean value of moments is (cf. (54)) 



i(l + 5)M(+) + i(l-5)M(-) 

[1 + (29-1 _ . _ S?r=o . ^tot(n ) Pj(^ 

(ELoi^tot(n)P„(£;)) 


(57) 


The last equation is valid in the first order on 
the only difference is the substitution 5 —>■ 5g. 
of 77-2 = 3 one has to use 77-2 = 3(1 + 5). 


5. one can performed exact summation and with the same accuracy 
In particular, the participation ratio is given by Eq. (50) but instead 


Inverse moments for the Ising model with the same values of parameters as above calculated numerically for zero 
momentum sector are plotted in Fig. for different values of q and compared with Eq. (57). Participation ratios for 
different values of a, A, and N are shown in Fig. In all cases in the bulk a good agreement with statistical model 
with corrections is found in the bulk of the spectrum. 
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FIG. 5. Participation ratio for the Ising model for indicated values of a and A and different N in zero translational momentum 
sector. From top to bottom N = 17, = 16, = 15, and N = 14. Small crosses are numerical results. White lines are the 

corrected statistical approximation. 


IV. CONCLUSION 


Exact determination of wave functions in many-body problems is rarely possible. For overwhelming majority of 
models numerical calculations remain the only way of getting information of wave functions. In certain models (called 
chaotic) the complexity of wave functions is so high that with a good precision they can be considered as random and 
the developing of statistical methods becomes a valuable alternative to numerics. This in turn enables meaningful 
formulations of questions about the foundations of quantum thermodynamics, which is a work in progress in modern 
physics. 

In the paper we discuss the construction of statistical models for eigenfunctions in a particular example of one¬ 
dimensional models of N interacting spin-^, namely the Ising model in transverse and longitudinal fields. The model 
is simple enough to demonstrate general phenomena without unnecessary complications. 

The investigation of the model is restricted to the chaotic regime when all coupling constants are of the same 
order. It is attested that for large number of spins eigenfunction coefficients in the bulk of the spectrum are well 
approximated by Gaussian functions with zero mean and variance determined analytically from the Hamiltonian. 
Such asymptotic results are supposed to give a good description of wave functions only in thermodynamic limit when 
the number of spins tends to infinity. For numbers of spins accessible in numerical calculations there exist small but 
noticeable deviations from asymptotic formulae and a large part of the paper is devoted to calculations of different 
types of corrections. 

The first type of corrections is common. It appears every time one calculates corrections to the central limit theorem 
by taking into account higher order moments. The best way to incorporate them is to use Gibbs-like formulae which 
requires numerical calculations to find necessary parameters. When the number of spins is large, contributions of 
higher order cumulants are small and one can incorporate these corrections by using simple series in the Hermite 
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polynomials. 

The second type of corrections appears for models with periodic boundary conditions and is related with the conser¬ 
vation of translational momentum and the parity transformation. Due to symmetry considerations some coefficients 
are real and others have to be complex. In addition, for certain values of translational momentum the symmetries 
impose that a subset of coefficients is zero. 

When all corrections are included, the resulting statistical expressions agree well with the results of numerical 
calculations of of wave function moments which clearly confirms that in the bulk wave functions are fully extended. 

All constructions used in the paper can be applied only in the bulk of the spectrum where eigen-energies are of the 
order of the square root of number of spins. The structure of eigenfunctions close to the ground state and to the highest 
energy state may and will be more complicated. In particular, the ground state eigenfunctions for practically all spin 
chains have non-trivial multifractal dimensions [35) . The existence or not a sharp transition between multifractal 
behaviour of boundary states and the extended states in the bulk and more careful determination of the nature of 
these states (glassy?) require additional investigations. The behaviour of wave functions for small or large values of 
coupling constants in the bulk is different from the one investigated so far. The existence in such cases of perturbation 
series for large but finite number of spins leads to multiple peaks in the both average spectral density [36j and in 
moments of wave functions. These questions will be discussed elsewhere. 


Appendix A: Calculations of first moments for the Ising model 


Hamiltonian of the Ising model in two fields ([^ can be rewritten in the form 

n = n, + n, (ai) 

where 

n, = -Y,K, h„ = «+i+a<, = (A2) 

n n 

The purpose of this Appendix is to calculate moments of the Hamiltonian iJ,k{n) = {n\'H'^\n) taken between basis 
states with definite projection of spins in each point. 

As Hzln) = En\n) with En = X{N — 2n) where n is the number of spins up, one has {n\'H’^\n) = E!^. Therefore the 
first two moments of the Ising model Hamiltonian are 

fii{n) = En, fi2{n) = N {1 + a'^) + El (A3) 

and the variance (i.e. the second cumulant) 

al = {n\'H‘^\n) - {n\'H\n)'^ = A^(l -I- a^) (A4) 

Hamiltonian T-Lx flips spins so simple counting gives 

{n\Hx\n)=0, (nlHlln) = N{1 + a"^), {n\Hl\n) =—6N (A5) 


The forth power of Hx is 

hl^hn ^+6 Y Y 

n ni^n2 ni<n2 ni^n2<nz 

+ 24 E ^ni hn2 ^714 

ni<n2<ra3<n4 

By inspection one can check that only the even powers of give contribution to (nlTf^ln) and 

{n\Hi\n) = 3A^2(1 -f + N{2Aa^ - 2 - 2a^) 

As it does not depend on n it coincides with 1/2-^Tr Hx 
The third power of the Hamiltonian is 


0/3 0/3 I o/2o/ \ U 'U 'U \'U 'U‘^ i U \ 'U 'U 'U i -TV^o/ i o/3 

Tl — f^x^z^x ~r f^x^z ' ^z^x ' ^z^x^z > ' ^z 

Due to Eq. (A5| one can find all contributions to {n\H^\fi) except HxHzHx- 


(A6) 


(A7) 


(A8) 
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One gets 


{n\n^\H) = {n\nl\n) + 2Er^{n\nl\n) + 3El{n\na=\n) + E^ + (A9) 

As has been discussed above, Hx acting on a product state |n) changes the number of spins up which can be 
symbolically written as follows 

T-Lx —> —1, 2/c times (AlO) 

n—^n 

T-Lx —>■ —a, (iV — n) times (All) 

n—¥n-\-l 

T-Lx — —ct, n times (A12) 

n—>-n —1 

Tix —>■ —1, (iV — n — A:) times (A13) 

n—>-n +2 

Tix —> —1, u — fc times (A14) 

n—^n—2 

where k is the number of groups of spins in the same direction. 

One gets 


T-izT-ix\n) = -2kEn\n) - [N -n - k)En+ 2 \n + 2) - (n - k)En- 2 \n - 2) 

— a{N — n)En+i\n + 1) — anEn-i\n — 1) (A15) 

Here \n) stays for a state with n spins up. 

Finally the action Tix will select the corresponding term and one obtains 

{n\TixTizTix\n) = 2kEn {N - n - k)En+2 + (n - k)En-2 + a'^[{N - n)En+i + nEn-i] 

= (N — n)En+2 + 'n-En-2 + — n)En+i + nEn-i] = [A^(l + a^) — 4 — 2a^]En (A16) 

Combining all terms together gives that for A^ > 4 

fisin) = {n\Ti^\n) = E^ + 2{N + k)En + {N - n - A:)A „+2 + {n - fc)A „_2 (A17) 

+ a^[2Er,N + {N - n)Er,+i + nA„_i - 6N] = E^ + [3(1 + a^)N - 4 - 2a^]Er, - 6Na^ 

Notice that the dependence on k disappears. 

Finally the third cumulant, k^ = fj ,3 — 3/i2/J.i + 2/xf, takes the following value 


fcg = -eNa"^ - 2En{a^+ 2) 
The next step is the calculation of the forth moment of Hamiltonian (lAl 


0/4 _ 0/4 I 0/00/ I o/zo/ o/ I o/io/z {Unj 07 ^ ^ U U U U U U _i_ U 

Tl — '^x ' '^x'^z ' '^x'^z'^x '^X'^Z'^X'^z ~r '^X'^z 

I 07 07 ^ _i_ 07 07 ^07 _i _07 07 07 07 _i _07 07 07 ^ i 07^072 i 07^07 07 _i_ 07^07 _i_ 07 ^ 

“r tvz^'r ~r ^z^'r^z < ^z^x^z^x ~r i^z^x^y ~r ~r i^yi^x^z ~r iLyiLx \ 


Lz ! l^X ' ^ ' ^Z' ^X' ^Z' ^X \ I Lz I l^x ' ^z ' '^z'^X' ' ^z ' ' ^^Z 

When calculating (n|'H‘^|n) many terms are known from the above expressions 

(fi|-H'‘|n) = (nj-H^jn) + 2 {n\Til\n) E^, + 3 {n\Til\n)El + A{n\Tix\n)El 
+ 2 {n\TixTizTix\n)En + + {f[\nlTi,Tix\n) + {n\TixTizTil\n) + {n\TixTilTix\n) 


(A18) 

(A19) 

(A20) 


Only three last quantities require separate calculation. _ 

Term {n\TixTilTix\n) can be performed as above and the result is the same as Eq. (A16) but with substitution 
Em —t 

{n\TixTilTix\n) = 2kE^ {N — n — fc)A ^_|_2 + (n — k)E '^_2 + a'^{N — n)E‘^_^_l + a'^nE^_i 
= -32k\^ + {N -S)El + 16Af[(iV -‘i)El+ dA^A^] (A21) 

Operator Tix flips one or two nearby spins, therefore only the following contribution are non-zero 

{n\TilTizTix\n) = -2a2(fj| ^ (XI'^" 2 ) 1^) 

n ni n 2 

-a^(n|(^<) (^<)h.(X«+i)K) (A22) 
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Consider the first term. Operator acting on |n) produces N — n states with n + 1 spins up and n states with 

n — 1 spins down (cf. ( [All I and (AlO)). Operator Hz multiplies them by En+i and En-i correspondingly. Other 
terms should combine to produce the initial state. In the end one gets 

{n\HlHzH4n) = -Aa\N - n)A„+i + nA„_i] (A23) 

— 2a^]p.kEn + {N — n — k)En +2 + {n — k)En- 2 ] = —2a^{3N — 8 )A„ 

Taking into account all terms we obtain that for N > 5 

k) = {n\K^\n) = E^ + El[6N{l + a^) - 16 - Sa^] + 8£;„(4 - aA^)^^ + 3N^{1 + 

- 32k\^ + N{2Aa^ -2-2a^ + 16X^ + dA^a^) (A24) 

The forth cumulant by definition is the following combination, ^4 = /i 4 — 4 /r 3 /ri — 3/^2 + 12 /r 2 /r 3 — . Simple 

calculations gives 


fc4 = 32a^En - 32kX^ + N{24:a'^ - 2 - 2a'^ + 16X^ + AX^o^) 


(A25) 


In the calculations we used instead of k its mean value of k in all states with fixed n, (k) = n{N — n)/{N — 1) [36] . 

All above calculations were performed without the selecting a particular translational momentum. If instead of basis 
set |n) one uses states which are eigenstates of the translational operator as in (30), moments may be different. We 
check that for lower order moments the difference is small and does not change noticeable the results. For simplicity 
we ignore such corrections. 
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